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COLISEUM  is  an  application  framework  that  integrates  plasma  propagation  schemes 
and  arbitrary  3D  surface  geometries.  Using  Particle-in-Cell  (PIC)  schemes  to  model  the 
plasma  propagation  high  fidelity  modeling  of  the  plasma  and  its  interactions  with  the  sur¬ 
faces  is  possible.  In  order  to  improve  the  computational  performance  of  the  Particle-in-Cell 
scheme  with  Direct  Simulation  Monte  Carlo  collision  modeling  (PIC-DSMC)  within  COL¬ 
ISEUM,  AQUILA,  acceleration  techniques  have  been  developed  that  significantly  decrease 
the  amount  of  CPU  time  needed  to  obtain  a  steady-state  solution.  These  techniques  have 
been  demonstrated  to  decrease  the  CPU  time  from  3  to  24  times  with  little  appreciable 
differences  in  the  global  particle  properties  and  number  densities.  This  work  investigates 
the  differences  in  the  local  plasma  properties  that  result  from  the  application  of  the  dif¬ 
ferent  acceleration  techniques.  Results  show  that  the  subcycling  acceleration  scheme  does 
accurately  capture  the  macroscopic  flow  properties  (such  as  particle  counts  and  species 
number  densities)  and  the  velocity  distributions  in  the  lower  density  regions  of  the  flow 
field.  However,  the  higher  density  regions  of  the  flow  field  (such  as  in  the  main  beam  of  the 
plasma  source)  show  significant  differences  that  are  believed  to  be  associated  with  the  sim¬ 
plifying  assumptions  used  in  the  original  collision  modeling  scheme  within  the  PIC-DSMC 
module  AQUILA. 


Nomenclature 


c  Velocity,  [m/s] 

Ci  Three-dimensional  velocity  space,  [m/s  x  m/s  x  m/s] 
cr  Relative  velocity  between  two  particles  [m/s] 

/  Velocity  distribution  function,  [-] 
k  Boltzmann  Constant,  [8.617339  x  10-5  eV/K] 

Mass  flow  rate,  [kg/s] 

Number  of  computational  particles 
n  Particle  number  density  [m~3] 

P  Probability  of  collision 

S  Surface  of  integration 

t  Time,  [s] 

Te  Electron  Temperature  [eV] 

V  Volume  of  computational  cell  [m3] 

Wp  Ratio  of  physical  particles  to  computational  particles 

x  Position,  [m] 
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Subscripts 

max  Maximum  value 
s  Surface  property 

Conventions 

DSMC  Direct  Simulation  Monte  Carlo  Collision  Modeling 

MCC  Monte  Carlo  Collision  Modeling 

NTC  No- Time-Counter  Collision  Model 

PIC  Particle-in-Cell  Modeling  Technique  for  Particles 

SCN  COLISEUM  cases  run  without  subcyling 

SCY  COLISEUM  cases  run  with  subcyling 

VDF  Velocity  Distribution  Function 

Symbols 

eo  Permittivity  of  free  space,  «  8.8542  x  10~12  F/m 

cf>  Electrostatic  potential  [V] 

pe  Electron  charge  density  [C/m3] 

pi  Ion  charge  density  [C/m3] 

a t  Total  collision  cross-section  [m2] 

r  Characteristic  time  [s] 


I.  Introduction 

The  Air  Force  Research  Laboratory  has  developed  an  application  framework  that  integrates  plasma  prop¬ 
agation  schemes  with  arbitrary  3D  surface  geometries1  in  order  to  investigate  the  interactions  between 
the  plasma  plume  from  electric  propulsion  thrusters  and  the  spacecraft.  A  hybrid  Particle-in-Cell  plasma 
model  within  COLISEUM,  called  AQUILA,2  is  the  basis  of  this  work.  The  COLISEUM  framework  is  used  to 
model  the  thrusters  in  space  environment  for  prediction  of  the  interactions  and  also  to  model  the  thrusters 
in  vacuum  chambers  in  order  to  validate  the  models  used  in  COLISEUM  against  experimental  data.  In 
order  to  improve  the  computational  performance  of  the  AQUILA  model,  acceleration  techniques  have  been 
developed  that  significantly  decrease  the  amount  of  CPU  time  needed  to  get  the  simulation  to  a  steady-state. 
These  schemes  have  been  demonstrated  to  decrease  the  CPU  time  by  up  to  24  times  with  little  appreciable 
differences  in  the  global  particle  properties. 

Two  previous  studies  have  been  performed  on  the  acceleration  schemes  within  COLISEUM.  Gibbons  et 
al.3  demonstrated  the  subcycling  technique  in  which  the  slower  moving  neutrals  are  propagated  at  a  larger 
time  step  than  the  faster  moving  ions.  Gibbons  et  al.4  demonstrated  the  use  of  the  subcycling  scheme  with  a 
scheme  that  decoupled  the  modeling  of  surface  interactions  from  the  plume  propagation.  It  utilizes  particle 
sources  from  the  surfaces  in  place  of  the  self-consistent  surface  interaction  modeling.  Both  of  these  schemes 
are  intended  to  provide  increased  convergence  rates  to  a  steady-state  solution.  The  desire  is  to  end  up  with 
the  same  final  plasma  distribution  with  or  without  the  acceleration  techniques.  The  present  study  intends 
to  provide  a  detailed  analysis  of  the  resulting  plasma  properties  and  any  differences  in  the  local  plasma 
properties  that  result  from  the  application  of  the  subcycling  acceleration  techniques. 

II.  Computational  Techniques 

This  section  provides  a  brief  overview  of  the  major  computational  techniques  within  COLISEUM  and 
AQUILA  obtained  from  the  User’s  Manuals. 

A.  COLISEUM  Framework 

COLISEUM  is  a  computational  framework  that  provides  a  tool  that  can  be  used  to  model  the  interaction 
of  plasmas  with  arbitrary  surfaces  in  three-dimensional  space.3  These  simulations  can  be  of  a  plasma  in  a 
contained  domain  or  in  an  open  domain.  This  allows  the  simulation  of  experiments  in  vacuum  chambers  as 
well  as  plasmas  in  the  low  density  space  environment.  The  primary  focus  of  COLISEUM  is  to  investigate 
the  erosion  associated  with  the  plasma  particles  impacting  surfaces,  known  as  sputtering,  as  well  as  the 
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re-deposition  of  this  material  on  other  surfaces.  COLISEUM  provides  the  integration  of  the  CAD  surface 
modeling,  the  plasma  propagation,  and  the  sputtering  of  material  within  one  consistent  framework. 

1.  Surface  Modeling 

The  surface  model  can  be  input  in  a  variety  of  standard  formats  including  ANSYS  and  ABAQUS  formats. 
In  addition,  surface  properties  are  also  specified  in  order  to  differentiate  the  various  materials  that  may 
compose  each  surface.  In  order  to  accurately  model  the  sputtering,  additional  material  information  must  be 
specified  that  describes  the  interaction  between  particles  that  may  be  impacting  surfaces  and  the  material 
that  the  surfaces  are  composed  of.  This  is  known  as  the  material  interaction  parameters  in  COLISEUM. 

The  sputtering  models  in  COLISEUM  are  based  on  standard  models  from  Roussel  et  al.5  and  Gardner 
et  al.,6  Kannenberg  et  al.,7  and  Yamamura  et  al.8  Coupling  the  sputtering  models  with  the  re-deposition 
process  has  been  included  in  order  to  account  for  how  the  re-deposited  material  may  itself  be  sputtered.  This 
allows  a  more  accurate  model  of  the  final  surface  deposition  characteristics.  Fife  et  al.1  have  also  developed 
an  iterative  scheme  to  model  resputtering  of  the  deposited  material. 


2.  Plasma  Modeling 

The  plasma  modeling  within  COLISEUM  has  two  major  components  to  it.  The  first  is  the  modeling  of  the 
source.  Sources  are  surfaces  within  the  geometry  from  which  particles  will  be  emitted.  To  model  sources,  a 
velocity  distribution  function,  VDF,  must  be  specified  throughout  the  surface  of  the  source.  In  general,  the 
VDF  is  a  function  of  space,  time,  and  velocity,  /  (x,  c,  t).  This  is  used  to  determine  the  mass  flow  rate,  m,  as 


rhs 


fs  (x,  c,  t)  dcdS 


(1) 


where  Cj  is  the  three-dimensional  velocity  space,  and  S  is  the  surface  of  the  source.  The  VDF  does  not  need 
a  specified  direction  since  COLISEUM  uses  the  surface  normals  from  the  geometry.  COLISEUM  provides 
surface  models  for  common  sources:  (1)  user  specified  flux  information  typically  from  experimental  data,  (2) 
user  specified  flux  and  velocity  information  typically  from  experimental  data,  and  (3)  a  shifted  Maxwellian 
distribution  with  the  velocity  shift  normal  to  the  surface. 

The  second  major  component  of  the  plasma  modeling  within  COLISEUM  is  the  plasma  simulation  itself. 
COLISEUM  was  developed  with  the  idea  of  supporting  any  number  of  plasma  modeling  schemes.  One  of  the 
original  models  is  the  prescribed  plume  model  which  simply  imports  previously  obtained  plasma  properties 
of  the  flow  field.  The  particle  fluxes  to  the  surfaces  are  calculated  by  mapping  the  solution  to  the  surface 
geometry.  The  other  original  model  is  a  ray  tracing  model  which  traces  the  particle  trajectory  from  the 
sources  without  accounting  for  the  electrostatic  potential  field  forces  or  particle  collisions.  The  particle 
fluxes  to  the  surfaces  are  then  determined. 

As  COLISEUM  has  matured,  more  sophisticated  plasma  modeling  modules  have  been  developed.  DRACO 
from  Virginia  Tech9  is  a  Cartesian  cell  based,  finite-element  Particle-in-Cell  Monte  Carlo  Collision  (PIC- 
MCC)  simulation.  AQUILA  from  MIT10  is  an  unstructured  tetrahedral  cell  based,  finite-element  PIC-DSMC 
simulation.  It  is  AQUILA  that  is  being  used  as  the  basis  of  this  investigation. 


B.  Particle  Propagation  Scheme 


To  propagate  the  particles,  two  separate  schemes  are  used.  The  neutrals  and  ions  are  treated  as  particles 
and  are  propagated  via  a  particle  tracking  technique.  The  electrons  are  treated  as  a  fluid.  The  electrons  are 
described  by  the  Boltzmann  relation 


ne 


\kTe) 


(2) 


The  time  integration  scheme  used  within  AQUILA  to  propagate  the  plasma  particles  is  the  standard  leap  frog 
scheme11  which  is  second  order  accurate  in  time.  The  electrostatic  forces  are  modeled  using  the  electrostatic 
potential  equation  with  the  inclusion  of  space-charge  effects10 


£o 


(3) 
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where  (f>  is  the  electrostatic  potential,  pe  is  the  electron  charge  density,  pi  is  the  ion  charge  density,  and  eo 
is  the  permittivity  of  free  space.  A  finite  element  formulation  is  used  to  solve  this  potential  equation  with  a 
Newton-method  type  scheme  to  handle  the  nonlinear  nature  of  the  resulting  equations. 


C.  Collision  Modeling 

The  collision  modeling  within  AQUILA12  is  based  on  the  No-Time-Counter,  NTC,  method  of  Bird.13  The 
probability  of  a  collision  between  two  particles  is  given  as 

p  =  WP  i?TCr)  A t  ^ 

where  Wp  is  the  ratio  of  physical  particles  to  computational  particles,  cft  is  the  total  collision  cross-section, 
cr  is  the  relative  speed  between  the  two  particles,  and  V  is  the  volume  of  the  computational  cell  containing 
particles.  Similarly,  the  maximum  probability  of  a  collision  is 


P  = 

'  max  — 


Wp  (°TCr)max  A  t 
V 


(5) 


The  NTC  scheme  samples  only  a  fraction  of  the  total  number  of  particle  pairs  in  the  computational  cell, 
and  adjusts  the  probability  of  collision  of  the  sampled  particle  pairs  accordingly.  Within  COLISEUM,  only 
Pmax Np Nq  particle  pairs  are  chosen  from  species  p  and  q.  Thus,  the  lower  the  maximum  probability  of 
collision,  the  fewer  collision  samples  are  taken.  The  resulting  collision  probability  for  a  sampled  collision 
pair  is  then 


P  = 


(7t  Cr 

(0"TCr)max 


(6) 


Notice  that  this  scheme  will  sample  the  appropriate  number  of  collision  pairs  only  when  an  accurate  maximum 
probability  has  been  determined.  Before  such  time,  too  few  collisions  pairs  will  be  sampled  resulting  in  the 
lower  probability  events  (but  not  insignificant  events)  being  under  represented.  Therefore,  this  scheme 
will  produce  accurate  collision  rates  only  after  a  sufficient  number  of  pairs  have  been  sampled  so  that  the 
maximum  probability  term  has  been  reasonably  determined.  Once  this  occurs  the  collision  calculations 
should  then  reasonably  capture  all  of  the  desired  collision  events. 

Within  AQUILA,  each  particle  species  (neutrals  and  ions  of  the  same  atom  are  considered  to  be  different 
species  in  addition  to  different  atoms)  can  have  its  own  physical  to  computational  particle  weighting.  This 
means  that  Wp  from  above  corresponds  to  the  larger  of  the  two  weightings  in  the  collision  pair.  Also,  while 
the  lower  weighted  pair,  Wq  in  the  collision  will  always  be  altered  by  the  collision,  when  one  occurs,  the 
probability  of  the  higher  weighted  particle  being  altered  by  the  collsion  is 


Pwp 


W i 
Wp 


(7) 


This  results  in  the  possibility  of  momentum  not  being  conserved  for  individual  collisions,  but  momentum 
will  be  conserved  on  average  with  a  sufficient  number  of  collisions. 

Finally,  for  all  collision  related  calculations,  a  simple  accept/reject  scheme  can  be  used  on  the  probabilities. 


D.  Subcycling  Scheme 

The  subcycling  scheme  within  COLISEUM  utilizes  the  fact  that  there  is  a  significant  difference  between  the 
collision  times  scales  and  the  ion  characteristic  time.  The  ion  characteristic  time  is  based  on  the  spatial 
resolution  of  the  computational  grid  which  is  being  used  to  model  the  electrostatic  potential  field.  In  order 
to  keep  the  electrostatic  forces  on  the  particles  varying  smoothly  as  they  travel  through  the  domain,  the 
ions  should  not  travel  more  than  a  third  of  a  cell  within  one  time  step.  For  the  simulation  to  be  modeled 
in  this  paper,  the  ion  velocity  is  20  km/s  and  the  neutral  velocity  is  200  m/s.  For  a  characteristic  length 
of  the  smallest  volume  of  0.01  m,  this  results  in  the  ion  characteristic  time  of  5  x  10-7  s  and  the  neutral 
characteristic  time  of  5  x  10-5  s.  This  is  a  factor  of  100  difference  between  the  two. 
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Using  the  following  relations  for  the  elastic  neutral-neutral  and  neutral-ion  collisions14 
the  collision  cross-sections 


to  characterize 


_el  2.117  x  10-18 

aXe-Xe  ~  o.24 

el  8.2807  x  10“16 

aXe-Xe+  —  “ 

and  for  the  charge-exchange  collisions  between  neutrals  and  ions15 

aXe-Xe+  =  1-1872  x  10_2°  [-23.3  log  (tv)  +  188.81] 
a  mean  time  to  collision  can  be  determined  as 


(8) 

(9) 


(10) 


1 

r  = - 

nacr 


(ii) 


Table  1  shows  the  characteristic  times  for  the  simulation  to  be  modeled  in  this  paper  and  a  maximum 
particle  number  density  of  1018  m-3.  Included  in  these  calculations  is  the  high  speed  neutrals  that  will  result 
from  previous  charge  exchange  collisions.  Their  number  density  will  be  shown  to  be  100  times  less  than  the 
maximum  particle  number  density.  Thus,  the  neutrals  in  the  collision  calculations  could  have  the  low  speed 
200  m/s  value  or  the  high  speed  20000  m/s  value.  Clearly  a  simulation  time  step  less  than  4.7  x  10-5  s  is 
needed  to  resolve  the  collision  time  scales. 


Table  1.  Collision  Characteristic  Times 


Collision  Type 

cr  [m/s] 

cr  [m2] 

T  [S] 

Xe-Xe  Elastic 

200 

5.94  x  10~19 

8.42  x  10"3 

Xe-Xe  Elastic 

20000 

1.97  x  10"19 

2.56  x  10"2 

Xe-Xe+  Elastic 

10000 

8.20  x  10"2° 

1.21  x  10’1 

Xe-Xe+  Elastic 

20000 

4.18  x  10"2° 

1.21  x  10"3 

Xe-Xe+  Charge  Exchange 

10000 

1.14  x  10"18 

8.81  x  10"3 

Xe-Xe+  Charge  Exchange 

20000 

1.05  x  10-18 

4.75  x  10"5 

Therefore,  there  is  only  one  physical  phenomena  that  requires  a  time  step  in  the  order  of  10_ '  s,  and 
that  is  capturing  the  electrostatic  forces  applied  to  the  ions. 

Computing  one  complete  computational  cycle  encompasses  the  following  steps 

1.  Subcycle  Fast  Particles 

(a)  Move  Fast  Particles 

(b)  Inject  Fast  Particles 

(c)  Update  Electrostatic  Fields 

2.  Move  Slow  Particles 

3.  Inject  Slow  Particles 

4.  Perform  Collisions 

By  only  moving  the  slow  particles  a  fraction  of  the  number  times  that  the  fast  particles  must  be  moved 
as  well  as  performing  the  collisions  on  the  coarse  time  step  a  significant  amount  of  computational  effort  is 
saved. 
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E.  Velocity  Distribution  Function  Probe 

In  order  to  determine  the  local  plasma  properties,  a  new  probe  was  introduced  into  the  AQUILA  probe 
architecture.  This  probe  samples  a  region  in  space  (currently  a  computational  cell)  and  stores  the  velocity 
of  every  particle  of  a  specified  type  that  resides  within  the  cell.  The  frequency  of  sampling  can  be  adjusted 
as  well  as  the  start  and  end  times  of  the  sampling.  Once  the  sampling  is  completed,  the  probe  sorts  the 
particles  into  bins  of  user  specified  sizes.  The  results  can  be  written  out  as  a  table  of  the  non-empty  bins,  or 
as  a  Tecplot  formatted  structured  grid  data  file.  Further  processing  is  possible  with  this  data  if  only  binning 
on  particle  speed  is  desired. 


III.  Results 

The  following  results  are  all  for  the  same  test  problem.  First,  solution  convergence  is  demonstrated  by 
using  a  fine  time  step  that  is  on  the  order  of  the  characteristic  time  step  of  the  electrostatic  forces.  It  is 
also  demonstrated  that  further  refinement  of  the  time  step  does  not  result  in  any  appreciable  change  in 
the  solution.  Second,  a  solution  is  presented  with  a  coarse  time  step  that  is  on  the  order  of  the  collision 
time  scale.  This  solution  is  compared  to  the  fine  time  step  solution.  Next,  a  solution  is  presented  using  the 
subcycling  scheme  discussed  above  with  the  ions  moving  at  the  fine  time  step  and  the  neutrals  moving  at 
the  coarse  time  step. 


A.  Test  Problem  Description 


The  test  problem  is  a  highly  simplified  geometry 
based  on  a  plasma  source  within  a  vacuum  chamber. 

Figure  1  shows  the  surface  meshes  associated  with 
the  test  problem.  The  plasma  source  is  a  small  cylin¬ 
der  with  the  cylinder  axis  aligned  with  the  z-axis. 

The  plasma  is  emitted  in  the  positive  z-direction 
from  that  particular  face  of  the  cylinder.  The  vac¬ 
uum  chamber  is  simplified  to  a  cylinder  with  the 
cylinder  axis  again  aligned  with  the  z-axis.  The 
plasma  source  is  firing  towards  one  end  of  the  cham¬ 
ber,  while  the  opposite  end  of  the  chamber  is  a  parti¬ 
cle  sink  such  that  any  particle  that  hits  that  surface 
leaves  the  computational  domain. 

The  chamber  has  a  diameter  of  1.5  m  and  a 
length  of  2  m.  The  plasma  source  has  a  diameter 
of  0.1  m  and  a  length  of  0.1  m.  The  distance  be¬ 
tween  the  plasma  source  face  and  the  chamber  face 
is  1.3  m 

The  plasma  source  is  composed  of  two  particle 
types,  low  speed  neutrals  and  high  speed  ions.  Both 
are  modeled  using  the  drifting  Maxwellian  source 
model  within  COLISEUM.16  The  neutral  drift  ve¬ 
locity  is  200  m/s  and  temperature  is  700  K.  The 

ion  drift  velocity  is  20  km/s  and  temperature  is  10  eV.  The  ratio  of  particle  weights  between  the  neutrals 
and  ions  is  300  :  1  and  the  physical  to  computational  particle  ratio  for  the  neutrals  is  6.00  x  1011.  Finally 
the  ion  mass  flow  rate  is  5.0  x  1CV7  kg/s,  and  the  neutral  mass  flow  rate  is  1.0  x  10~7  kg/s. 

The  electrostatic  potential  is  modeled  using  the  quasi-neutral  model  within  AQUILA12  instead  of  the 
non-equilibrium  model  mentioned  previously.  This  applies  the  quasi-neutrality  assumption  and  inverts  Boltz¬ 
mann’s  equation  to  obtain  an  expression  for  the  electrostatic  potential 


Figure  1.  Simple  Chamber  Geometry 
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where  the  electron  temperature,  Te,  is  set  to  2  eV.  The  reference  electron  number  density  and  potential  is 
specified  to  be  at  a  potential  of  0  V  just  in  front  of  the  thruster  face. 
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In  order  to  examine  the  similarities  of  the  local  properties  of  the  plasma  between  the  three  cases,  the 
velocity  distribution  for  the  neutrals  and  ions  were  obtained  0.1  m  in  front  of  the  plasma  source  as  well  as 
0.28  m  above  the  thruster  face.  The  first  sampling  will  examine  the  plume  modeling  capabilities,  while  the 
second  sampling  will  examine  the  capabilities  to  model  the  plasma  outside  of  the  main  plasma  beam. 

The  computer  that  performed  these  simulations  for  the  timing  results  is  a  dual  processor  AMD  Opteron 
242  system  with  2  GB  of  RAM  with  an  additional  2  GB  of  swap  space.  Use  of  the  machine  was  minimized 
while  the  cases  were  running,  and  all  cases  resided  in  physical  memory,  so  the  swap  space  was  not  utilized 
except  to  move  other  non-essential  applications  out  of  RAM  at  the  start  of  the  simulation. 

In  order  to  quickly  distinguish  between  the  various  cases  to  be  run,  the  acronym  SCN  will  designate  cases 
without  the  use  of  subcycling  and  SCY  will  designate  the  cases  with  subcycling. 

B.  Solution  Convergence  Demonstration 

With  the  baseline  fine  time  step  for  this  simulation  established  as  2.5  x  10”'  s,  three  cases  were  run  to 
demonstrate  the  convergence  of  the  solution  at  this  time  step.  One  case  was  at  twice  the  baseline  time 
step,  5.0  x  10”7  s,  and  one  at  half  the  baseline  time  step,  1.25  x  10”'  s.  Each  case  was  computed  to  a 
final  computational  time  of  0.25  s.  Table  2  shows  the  collision  rates  for  the  three  cases.  There  is  very  little 
difference  between  all  of  the  collision  rates  for  the  three  cases.  The  differences  between  the  cases  is  most 
certainly  due  to  the  statistical  scatter  associated  with  these  types  of  schemes.  Notice  that  halving  or  doubling 
the  step  size  approximately  doubled  or  halved  the  amount  of  time  required  to  calculate  the  solution. 

Table  2.  Collision  Rates  for  Solution  Convergence  Cases 


Scheme 

Compute 
Time  [hr] 

Time 

Step  [s] 

Total 

Collisions  [#/s] 

Xe-Xe+  Charge 
Exchange  [#/s] 

Xe-Xe 
Elastic  [#/s] 

Xe-Xe+ 
Elastic  [#/s] 

SCN 

107.2 

1.25  x  10”7 

5.301  x  107 

5.046  x  107 

2.519  x  105 

2.293  x  106 

SCN 

58.9 

2.5  x  10”7 

5.325  x  107 

5.069  x  107 

2.560  x  105 

2.303  x  106 

SCN 

30.9 

5.0  x  10”7 

5.310  x  107 

5.056  x  107 

2.541  x  105 

2.287  x  106 

Figure  2  shows  the  evolution  of  the  total  number  of  neutrals  and  ions  as  well  as  the  total  number  of 
particles.  These  counts  are  nearly  identical  with  any  differences  with  the  statistical  scatter. 


0  0.05  0.1  0.15  0.2  0.25 

time  (s) 

(c)  SCN  At  =  5.0  x  10”7  s 


Figure  2.  Global  Particle  Counts  for  Solution  Convergence  Cases 


Figure  3  shows  a  contour  plot  of  the  final  number  density  distribution  for  the  neutrals  for  the  three  cases. 
The  plane  that  is  shown  is  the  x-z  plane  through  center  of  the  plasma  source.  All  figures  shown  are  for  the 
same  plane  at  the  same  location.  Figure  4  shows  a  contour  plot  of  the  final  number  density  distribution  for 
the  ions  for  the  three  cases,  and  Figure  5  shows  a  contour  plot  of  the  final  electrostatic  potential  for  the  three 
cases.  All  of  these  figures  also  show  only  very  minor  differences  between  the  three  cases.  The  outer  wings 
of  the  plume  are  captured  in  all  three  cases  as  can  be  seen  in  the  ion  number  density  and  the  electrostatic 
potential.  In  addition,  the  time  step  is  clearly  sufficient  to  accurately  capture  the  ion  trajectories.  As  can 
be  seen,  the  ions  are  occupying  the  region  behind  the  plasma  source  which  can  only  occur  if  the  time  step  is 
sufficiently  small  to  allow  the  electrostatic  forces  to  curve  the  trajectories  of  the  charge-exchange  ions  away 
from  the  plasma  source.  The  very  low  ion  density  immediately  behind  the  ion  source  is  due  to  the  ions 
colliding  with  that  surface  and  reflecting  back  as  accommodated  neutrals. 
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Figure  3.  Final  Neutral  Number  Density  for  Solution  Convergence  Cases 
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(a)  SCN  At  =  1.25  x  10~7  s 


(b)  SCN  At  =  2.5  x  10“7  s 


(c)  SCN  At  =  5.0  x  10"7  s 


Figure  4.  Final  Ion  Number  Density  for  Solution  Convergence  Cases 


(a)  SCN  At  =  1.25  x  10“7  s 


(b)  SCN  At  =  2.5  x  10-7  s 


(c)  SCN  At  =  5.0  x  10-7  s 


Figure  5.  Final  Electrostatic  Potential  for  Solution  Convergence  Cases 
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Next,  the  local  properties  of  the  plasma  will  be  compared  between  the  three  cases.  The  sampling  is 
performed  using  the  VDF  probe.  Samples  are  obtained  every  2.5  x  10~5  s,  for  this  case  every  100  time  steps. 
Sampling  is  started  after  0.4  s  have  elapsed,  when  the  global  particle  count  has  reached  a  nearly  steady  state. 
The  neutral  sampling  within  the  plume  is  shown  in  Figure  6  for  both  the  velocity  distribution  as  well  as 
the  speed  distribution.  Thehe  velocity  distribution  is  showing  the  outer  edge  of  the  velocity  space  domain 
that  is  populated  with  particles  and  the  surfaces  are  colored  by  z-velocity  since  that  is  the  primary  velocity 
direction.  There  are  a  total  of  approximately  6100  particles  being  sampled  in  each  case.  Figure  7  shows  the 
ion  sampling  within  the  plume,  where  a  total  of  approximately  670,000  particles  are  being  sampled  in  each 
case.  Once  again,  there  is  very  little  differences  between  the  three  cases.  A  bimodal  distribution  is  seen  in 
the  ion  distribution  with  the  lower  speed  ions  being  the  result  of  the  charge-exchange  collisions  between  the 
neutrals  and  the  ions.  The  corresponding  high  speed  neutrals  can  be  seen  in  the  neutral  distributions,  but 
the  clarity  is  obscured  because  of  the  relative  low  occurrence  of  these  particles  due  to  the  computational 
particle  weighting  used  for  the  neutrals. 


(a)  SCN  Velocity  At  =  1.25  x  10  7  s 


speed  (m/s) 


(d)  SCN  Speed  At  =  1.25  x  10"7  s 


(b)  SCN  Velocity  At  =  2.5  x  10  7  s 


speed  (m/s) 


(e)  SCN  Speed  At  =  2.5  x  10  7  s 


(c)  SCN  Velocity  At  =  5.0  x  10  7  s 


speed  (m/s) 


(f)  SCN  Speed  At  =  5.0  x  1CT7  s 


Figure  6.  Neutral  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Solution  Convergence  Cases 


The  neutral  sampling  outside  the  plume  is  shown  in  Figure  8  for  both  the  velocity  distribution  as  well  as 
the  speed  distribution,  with  the  same  sampling  configuration  parameters  used  for  the  previous  samplings. 
Figure  9  shows  the  ion  sampling  outside  the  plume.  In  this  case  there  is  very  little  difference  between  the 
three  cases,  but  some  minor  differences  do  occur.  The  majority  of  the  neutrals  that  are  being  sampled  in 
this  region  are  due  to  the  ions  accomodating  and  reflecting  off  of  the  walls.  That  is  why  the  most  probable 
velocity  is  so  low.  However,  the  5.0  x  10-'  s  case  does  show  some  added  noise  in  the  higher  velocities.  Notice 
that  at  this  time  step,  we  are  very  close  to  the  ion  characteristic  time  step,  so  the  increase  in  the  time  step 
size  might  have  altered  the  ion  trajectories  enough  to  alter  the  neutral  distribution  in  the  higher  speeds. 
There  still  exists  a  bimodal  distribution  in  the  ion  velocity,  but  the  number  of  low  speed  ions  is  significantly 
less  than  what  were  in  front  of  the  thruster,  as  expected.  Also,  the  most  probable  speed  of  the  ions  has 
dropped  from  the  20  km/s  that  it  was  in  front  of  the  thruster  to  around  4  km/s.  This  is  because  the  ions 
that  make  it  to  this  region  of  the  flow  are  the  low  speed  charge  exchange  ions  that  have  been  accelerated 
through  the  electrostatic  potential  held  into  this  region. 
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(a)  SCN  Velocity  At  =  1.25  x  10  7  s 


(b)  SCN  Velocity  At  =  2.5  x  10  7  s 


(c)  SCN  Velocity  At  =  5.0  x  10  7  s 


speed  (m/s) 

(d)  SCN  Speed  At  =  1.25  x  10"7  s 


speed  (m/s) 

(e)  SCN  Speed  At  =  2.5  x  10— 7  s 


speed  (m/s) 

(f)  SCN  Speed  At  =  5.0  x  1CT7  s 


Figure  7.  Ion  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Solution  Convergence  Cases 


Figure  8.  Neutral  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Solution  Convergence 
Cases 
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(d)  SCN  Speed  At  =  1.25  x  lO"7  s 


Figure  9.  Ion  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Solution  Convergence  Cases 


C.  Coarse  Time  Step  Solution 

Now  the  coarse  time  step  case  can  be  compared  to  the  fine  time  step.  It  is  understood  that  these  results  will 
have  significant  errors  associated  with  the  time  step  being  too  large  to  capture  the  correct  ion  trajectories. 
This  case  is  still  instructive  in  observing  what  features  are  not  being  captured  by  the  coarse  time  step  and 
for  future  discussions  of  features  in  the  subcycling  scheme  that  can  be  attributed  to  the  coarse  time  step 
used  for  the  neutral  time  scale. 

For  this  case  all  particles  and  physical  processes  are  propagated  at  a  time  step  of  2.5  x  10-5  s.  From  the 
previous  discussions,  it  is  known  that  this  time  step  is  larger  than  the  characteristic  time  associated  with 
the  electrostatic  forces,  5.0  x  10~7  s,  and  is  very  close  to  the  smallest  characteristic  time  associated  with  the 
collision  modeling,  4.7  x  10~5  for  most  of  the  charge  exchange  collision  events  that  occur.  Table  3  shows  the 
resulting  collision  rates  along  with  the  fine  time  step  collision  rates. 

Table  3.  Collision  Rates  for  Coarse  Time  Step 

Compute  Time  Total  Xe-Xe+  Charge  Xe-Xe  Xe-Xe + 

Scheme  Time  [hr]  Step  [s]  Collisions  [#/s]  Exchange  [#/s]  Elastic  [#/s]  Elastic  [#/s] 
SCN  58.9  2.5  x  10~7  5.325  x  107  5.069  x  107  2.560  x  105  2.303  x  106 

SCN  1.73  2.5  x  10-5  4.945  x  107  4.726  x  107  2.572  x  105  1.936  x  106 


The  first  thing  to  notice  from  this  table  is  that  taking  100  times  fewer  time  steps  results  in  a  significant 
decrease  in  compute  time,  by  a  factor  of  around  34.  Unfortunately,  only  the  Xe-Xe  elastic  collision  rate  is 
the  same  between  the  coarse  and  fine  time  steps.  The  elastic  and  charge  exchange  collisions  associated  with 
the  Xe-Xe+  pairs  differ  between  the  two  cases.  Therefore  by  not  adequately  resolving  the  ion  trajectory 
there  has  been  a  decrease  in  the  collision  rate  associated  with  the  ions.  This  could  be  caused  by  the  ions 
traveling  entirely  through  the  high  density  region  in  front  of  the  plasma  source  (where  collisions  are  most 
likely)  before  the  ions  can  participate  in  a  significant  number  of  collision  events.  Taking  the  nominal  ion 
velocity  of  20  km/s  and  the  time  step  of  2.5  x  10-5  s  yields  a  distance  traveled  by  an  ion  of  0.5  m.  This  is  a 
significant  distance  from  the  plasma  source  and  is  certainly  outside  of  the  high  density  region  found  in  the 


11  of  23 


American  Institute  of  Aeronautics  and  Astronautics  Paper  2006-3248 
Distribution  A:  Approved  for  public  release;  distribution  unlimited 


fine  time  step  cases  from  Figures  3  and  4. 

While  significant  differences  are  seen  in  the  collision  rates,  Figure  10  shows  that  the  evolution  of  the  total 
number  of  neutrals  and  ions  as  well  as  the  total  number  of  particles  is  very  similar.  Thus,  the  differences 
between  the  two  cases  must  be  for  only  a  small,  but  significant,  fraction  of  the  total  particles. 


(a)  SCN  At  =  2.5  x  10"7  s  (b)  SCN  At  =  2.5  x  10"5  s 

Figure  10.  Global  Particle  Counts  for  Coarse  Time  Step 


Figure  11  shows  the  final  number  density  for  the  neutrals  for  the  coarse  and  fine  time  step  cases.  Even 
for  this  case  the  neutral  number  density  is  fairly  similar.  It  appears  that  the  coarse  time  step  does  not  have 
a  significant  effect  on  the  overall  neutral  number  density.  This  does  not  mean,  however,  that  there  is  no 
effect.  Since  the  collision  rate  is  different,  there  is  likely  some  difference  in  the  velocity  distribution  function 
throughout  the  computational  domain. 
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(a)  SCN  At  =  2.5  x  10“7  s 


(b)  SCN  At  =  2.5  x  lCT5  s 


Figure  11.  Final  Neutral  Number  Density  for  Coarse  Time  Step 


Figure  12  shows  the  final  number  density  for  the  ions  for  the  coarse  and  fine  time  step  cases.  This  shows 
a  significant  difference  between  the  two  cases.  While  the  main  beam  region  seems  similar,  the  outer  wings  of 
the  plume  are  certainly  not  captured  as  well  in  the  coarse  time  step  case.  Also,  with  the  time  step  so  large, 
the  ions  cannot  make  the  curved  trajectory  to  collide  with  the  back  of  the  ion  source,  which  results  in  the 
difference  between  the  two  cases  in  that  region. 

Figure  13  shows  the  final  electrostatic  potential  for  the  coarse  and  fine  time  step  cases.  This  again  shows 
less  accurate  resolution  of  the  outer  wings  of  the  plume  region.  Also  an  increased  potential  region  behind 
the  plasma  source  exists  since  the  ions  are  not  colliding  with  the  back  of  the  plasma  source  and  becoming 
neutralized. 

Next,  the  local  properties  of  the  plasma  will  be  compared  between  the  three  cases.  The  neutral  sampling 
within  in  the  plume  is  shown  in  Figure  14  for  both  the  velocity  distribution  as  well  as  the  speed  distribution. 
The  most  significant  difference  here  is  that  there  is  a  significantly  larger  band  of  high  velocity  neutrals, 
around  the  20-25  km/s  range,  with  a  noticeable  decrease  in  population  in  the  10  krn/s  region  for  the  coarse 
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Figure  12.  Final  Ion  Number  Density  for  Coarse  Time  Step 
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Figure  13.  Final  Electrostatic  Potential  for  Coarse  Time  Step 
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time  step.  The  20  km/s  neutrals  are  results  of  the  charge  exchange  collisions  with  the  high  speed  ions.  The 
10  km/s  neutrals  are  either  from  secondary  collisions  between  the  high  speed  neutral  and  the  other  particles 
(such  as  elastic  collisions  between  ions  and  neutrals  and  two  neutrals)  or  from  ellastic  collisions  between  the 
beam  ions  and  the  beam  neutrals.  One  point  to  note  is  that  with  the  large  time  step,  the  high  speed  neutrals 
are  able  to  travel  entirely  through  the  high  density  region  before  any  other  collision  event  might  occur. 
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Figure  14.  Neutral  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Coarse  Time  Step 

Figure  15  shows  the  ion  sampling  within  the  plume.  Again  a  difference  can  be  seen  between  the  coarse 
and  fine  time  step  cases.  First,  there  is  a  noticable  increase  in  the  number  of  ions  in  between  the  two  modes 
in  the  distribution.  This  velocity  range,  around  10  km/s,  corresponds  to  the  decrease  in  the  distribution  of 
neutrals  mentioned  above.  It  appears  that  with  the  coarse  time  step  there  is  not  enough  of  the  secondary 
collisions  to  decrease  the  population  of  the  higher  energy  neutral  as  well  as  the  moderate  energy  ions.  Notice 
that  these  ions  are  mainly  produced  by  the  charge  exchange  collisions  from  a  previous  time  step  since  the 
plasma  source  is  producing  a  Maxwellian  distribution  of  ions  with  the  peak  of  the  distribution  at  20  km/s. 
Thus  it  appears  that  the  coarse  time  step  is  moving  the  particles  out  of  the  high  density  region  too  quickly. 

The  neutral  sampling  outside  the  plume  is  shown  in  Figure  16  for  both  the  velocity  distribution  as  well 
as  the  speed  distribution.  While  the  two  distributions  look  similar,  there  is  a  more  concentrated  collection 
of  particles  in  the  8  km/s  region  for  the  coarse  time  step.  This  is  hard  to  distinguish  from  the  statistical 
scatter  in  the  data,  but  this  clustering  does  occur  over  8  consecutive  velocity  bins,  so  this  is  likely  more  than 
a  statistical  artifact. 

Figure  17  shows  the  ion  sampling  outside  the  plume.  In  this  case  there  are  significant  differences  between 
the  coarse  and  fine  time  step  cases.  The  most  notable  is  that  the  coarse  time  step  has  a  bimodal  distribution 
with  the  second  peak  rather  wide  and  centered  around  10  km/s.  This  peak  does  not  appear  on  the  fine  time 
step  and  must  be  a  result  of  the  coarse  time  step.  Again,  the  speed  is  associated  with  secondary  collisions, 
and  if  the  fast  moving  particles,  which  are  the  only  ones  that  create  these  particles,  travel  through  the  high 
density  region  before  another  collision  event  is  performed,  then  there  would  be  left  over  medium  speed  ions. 
Also,  the  width  of  the  main  peak,  that  is  center  at  4  km/s  is  much  larger  for  the  coarse  time  step  than  for 
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Figure  15. 


Ion  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Coarse  Time  Step 
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Figure  16.  Neutral  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Coarse  Time  Step 


16  of  23 


American  Institute  of  Aeronautics  and  Astronautics  Paper  2006-3248 
Distribution  A:  Approved  for  public  release;  distribution  unlimited 


the  fine  time  step.  There  is  also  a  corresponding  increase  in  the  speed  distribution  function  value  at  the 
peak.  Again  since  this  peak  is  most  likely  a  result  of  multiple  collisions  occuring  before  the  particle  leaves 
the  high  density  region,  the  coarser  time  step  can  again  account  for  this  difference. 
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Figure  17.  Ion  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Coarse  Time  Step 


It  appears  that  there  is  a  significant  coupling  between  the  particle  time  step  and  the  ability  to  capture 
all  of  the  secondary  collisions  that  are  occurring.  Even  though  the  time  step  was  fine  enough  for  the  collision 
modeling  characteristic  time  step,  it  appears  that  the  rapid  particle  density  variation  in  front  of  the  plasma 
source  is  causing  the  collision  characteristic  time  step  calculation  to  be  too  large.  Also,  notice  that  while 
some  of  these  effects  are  noticeable  in  the  gross  perspective  of  the  flow  field,  the  significant  differences  are 
only  apparent  with  the  observation  of  the  velocity  distribution  functions. 


D.  Subcycling  Solution 

Now  that  is  apparent  the  the  coarse  time  step  does  not  capture  a  number  of  significant  flow  features,  an 
analysis  of  the  improvements  associated  with  the  subcycling  algorithm  can  be  performed.  The  subcycling 
algorithm  uses  the  coarse  time  step,  2.5  x  10-5  s,  for  the  slow  particle  time  step  (i.e.,  slow  neutrals)  and  the 
fine  time  step,  2.5  x  1CD7  s,  for  the  fast  particle  time  step  (i.e.,  for  the  fast  ions  and  neutrals).  Table  4  shows 
the  resulting  collision  rates  along  with  the  fine  and  coarse  time  step  collision  rates  for  comparison. 

The  first  thing  to  notice  from  this  figure  is  that  the  subcycling  scheme  results  in  a  decrease  of  computa¬ 
tional  time  from  58.9  hr  to  13.7  hr  compared  to  the  fine  time  step  case.  Unfortunately,  the  collision  rates 
are  nearly  identical  to  the  coarse  time  step  rates  and  are  significantly  lower  than  the  fine  time  step  results 
(with  the  same  exception  of  the  Xe-Xe  elastic  collision  rate). 

While  significant  differences  are  seen  in  the  collision  rates,  Figure  18  shows  the  evolution  of  the  total 
number  of  neutrals  and  ions  as  well  as  the  total  number  of  particles  is  very  similar.  Thus,  the  differences 
between  the  subcycling  case  and  the  fine  time  step  case  again  must  be  for  only  a  small,  but  significant, 
fraction  of  the  total  particles. 
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Table  4.  Collision  Rates  for  Subcycling  Case 
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Figure  18.  Global  Particle  Counts  for  Subcycling  Case 


Figure  19  shows  the  final  number  density  for  the  neutrals  for  the  subcycling  case  as  well  as  the  fine  and 
coarse  time  step  cases.  This  case  again  shows  that  there  is  little  difference  between  the  neutral  number  den¬ 
sities  for  the  subcycling  case.  Thus,  the  overall  neutral  distribution  is  fairly  insensitive  to  the  computational 
time  step.  However,  as  was  the  case  for  the  coarse  time  step,  it  is  expected  that  there  will  be  differences 
that  are  not  observable  in  this  perspective. 


Figure  19.  Final  Neutral  Number  Density  for  Subcycling  Case 


Figure  20  shows  the  final  number  density  for  the  ions  for  the  subcycling  case  as  well  as  the  fine  and 
coarse  time  step  cases.  Unlike  the  coarse  time  step  case,  the  ion  number  density  is  very  close  to  the  fine  time 
step  case.  The  outer  wings  of  the  plume  are  captured,  and  the  ion  neutralization  at  the  back  of  the  plasma 
source  is  also  captured.  Therefore,  the  subcycling  is  drastically  improving  the  capabilities  of  capturing  the 
ion  spatial  distribution. 

Figure  21  shows  the  final  electrostatic  potential  for  the  subcycling  case  as  well  as  the  fine  and  coarse  time 
step  cases.  This  again  shows  that  the  subcycling  case  and  the  fine  time  step  cases  are  quite  similar.  This  is 
to  be  expected  since  the  electrostatic  potential  is  directly  related  to  the  ion  distribution.  Even  the  potential 
drop  behind  the  thruster  seen  in  the  fine  time  step  case  is  captured  in  the  subcycling  case. 

While  the  collision  rates  are  different,  the  subcycling  case  has  so  far  improved  the  ion  number  density 
distribution  compared  to  the  coarse  time  step  case  and  has  shown  no  difference  in  the  neutral  number  density 
distribution. 
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Figure  20.  Final  Ion  Number  Density  for  Subcycling  Case 


Figure  21.  Final  Electrostatic  Potential  for  Subcycling  Case 
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Next,  the  local  properties  of  the  plasma  will  be  compared  between  the  subcycling  case  and  the  fine 
and  coarse  solution.  The  neutral  sampling  within  the  plume  is  shown  in  Figure  22  for  both  the  velocity 
distribution  as  well  as  the  speed  distribution.  Unfortunately,  the  subcycling  cases  looks  much  more  similar 
to  the  coarse  time  step  case  and  has  significant  differences  with  the  fine  time  step  case.  The  same  arguments 
about  the  coarse  time  step  differences  also  seem  to  apply  here.  While  the  neutrals  are  propagating  at  the 
fine  time  step,  there  is  still  no  mechanism  to  get  these  neutrals  to  participate  in  collision  events  while  they 
reside  in  the  high  density  regions  since  collisions  are  only  computed  at  the  coarse  time  scale. 
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Figure  22.  Neutral  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Subcycling  Case 


Figure  23  shows  the  ion  sampling  within  the  plume.  This  subcycling  case  also  shows  significant  differences 
between  the  subcycling  and  fine  time  steps  and  is  very  similar  to  the  coarse  time  step.  It  appears  that  this, 
too  might  be  attributed  to  the  secondary  collisions  discussed  previously. 

The  neutral  sampling  outside  the  plume  is  shown  in  Figure  24  for  both  the  velocity  distribution  as  well 
as  the  speed  distribution.  While  it  is  not  certain  that  the  8  km/s  region  in  the  coarse  time  step  is  caused  by 
statistical  scatter,  it  is  worth  noting  that  the  subcycling  case  does  not  demonstrate  this  feature.  Otherwise, 
the  subcycling  case  looks  very  similar  to  the  fine  time  step  case. 

Figure  25  shows  the  ion  sampling  outside  the  plume.  This  case  shows  drastic  improvements  from  the 
coarse  time  step.  The  subcycling  case  does  not  have  the  secondary  peak  in  the  10  km/s  range  and  has 
a  similarly  narrow  speed  range  around  the  most  probable  speed.  It  is  apparent  that  the  subcycling  does 
significantly  improve  the  particle  modeling  outside  the  high  density  plume  region.  This  is  most  likely  due 
to  the  fact  that  the  trajectory  of  the  high  speed  ions  is  significantly  improved  with  the  fine  time  step,  and 
thus  these  high  speed  ions  are  more  effected  by  the  electrostatic  potential  field. 

IV.  Conclusions 

The  subcycling  acceleration  scheme  within  the  AQUILA  plasma  modeling  module  of  COLISEUM  was 
investigated  to  determine  how  effective  it  is  in  capturing  the  local  plasma  properties.  First,  the  simulation 
was  demonstrated  to  be  capable  of  converging  to  a  solution  for  a  sufficiently  fine  time  step.  This  solution 
was  the  used  to  compare  the  performance  of  the  simulation  at  a  much  coarser  time  step.  This  showed  several 
deficiencies  in  the  coarse  time  step  solution.  These  were  mainly  focused  on  the  fact  that  the  high  speed 
particles  are  leaving  the  high  density  region  where  multiple  collisions  are  expected  to  occur  after  one  or  two 
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(a)  SCN  Velocity  At  =  2.5  x  10  '  s 


(b)  SCN  Velocity  At  =  2.5  x  10  5  s 


(c)  SCY  Velocity  At  =  2.5  x  10  5  s 


Figure  23.  Ion  Velocity  and  Speed  Distributions  0.1  m  in  Front  of  Thruster  for  Subcycling  Case 
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Figure  24. 


Neutral  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Subcycling  Case 
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Figure  25.  Ion  Velocity  and  Speed  Distributions  0.28  m  Above  Thruster  Face  for  Subcycling  Case 


time  steps.  In  addition  the  coarse  time  step  resulted  in  the  ion  trajectories  being  significantly  off  which 
resulted  in  large  differences  in  the  ion  density  distribution  compared  to  the  fine  time  step  case. 

The  subcycling  case  improved  the  modeling  of  the  number  densities  of  the  ions  and  neutrals  compared 
to  the  coarse  time  step  simulation  with  a  compute  time  speedup  of  a  factor  of  3.4.  It  also  significantly 
improved  the  modeling  of  the  lower  density  region  outside  the  main  plume.  In  this  region  the  velocity 
distribution  functions  for  the  neutrals  and  ions  were  both  very  similar  to  the  fine  time  step  case.  This  shows 
that  improving  the  modeling  of  the  electrostatic  forces,  via  finer  time  steps  to  propagate  the  ions  and  the 
force  calculations,  does  improve  the  modeling  of  the  plasma.  For  the  higher  density  region  of  the  main  beam 
of  the  plasma,  the  subcycling  showed  no  improvement  compared  to  the  coarse  time  step  solution.  While 
this  appears  to  be  rooted  in  the  fact  that  a  high  speed  particle  can  travel  through  the  entire  high  density 
plume  region  in  one  time  step  and  thus  not  participate  in  the  additional  collisions  that  appear  to  produce 
the  different  features  of  the  fine  time  step  solution.  An  alternative  explanation  to  this  behavior  could  be 
that  the  collision  model  is  not  providing  consistent  results  at  the  different  time  scales.  Suppose  the  collision 
model  was  selecting  too  many  collisions  to  occur  at  the  finer  time  step  because  the  selection  criteria  was 
slightly  off.  This  would  result  in  a  non-physical  increase  in  additional  collisions  that  the  high  speed  particles 
would  be  participating  in.  This  would  mean  that  the  fine  time  step  solution  was  not  correct. 

It  is  apparent  from  this  study  that  the  collision  model  used  with  AQUILA  needs  further  validation. 
Although  AQUILA  has  been  successful  in  previous  studies  that  examined  the  aggregate  properties,2, 3’ 10, 12 
a  more  detailed  examination  of  the  collision  algorithm  is  needed.  The  combination  of  using  this  algorithm 
in  which  variable  species  weighting  is  used  with  the  subcycling  technique,  i.e.  separate  time  steps  for  high 
speed  and  low  speed  particles,  requires  more  attention.  For  this  type  of  rarefield  flow  in  which  collisions 
result  in  orders  of  magnitude  changes  in  velocity,  it  may  be  necessary  to  use  a  different  collision  algorithm 
and/or  conserve  momentum  and  energy  for  each  collision.  A  future  study  of  the  effects  of  various  collision 
algorithms  on  the  local  velocity  distribution  function  should  be  performed  to  assess  their  ability  to  handle 
the  variable  species  weighting  that  is  used  within  AQUILA. 
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